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Free energy analysis for a system of interacting particles arranged in Bravais lattice 
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We propose a method of free energy calculation for a system of interacting particles arranged in 
a Bravais lattice. It will be shown how to treat divergences for infinite unbounded systems with 
“catastrophic” potentials like Coulomb and compare two of them. Besides, we show that current 
method may be used not only for essentially classical systems, but for some quantum as well. 
Two systems are considered: grains in dusty plasma and electrons on the liquid-helium surface. 
For dusty-plasma parameters of grain’s lattice are calculated by numerical solution of obtained 
equations. Electrons on the liquid-helium surface are analyzed to get localization distance. Besides, 
conditions for existence of electron lattice are found. 
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INTRODUCTION 


There are many soft-matter systems, such as grains 
in dusty plasma, colloids in various solvents, surfactant 
solutions, etc., that exhibit self-organization and rear¬ 
rangement in crystalline structures under certain condi¬ 
tions. Often their inter-particle potential is long-range 
and similar to the Coulomb one. Usually, this kind of 
potentials cause divergences (that’s why they are often 
called “catastrophic”) and complicates consideration [ll- 
1^. On the other hand, mentioned systems are very inter¬ 
esting due to their application to the studies of a varieW 
of peculiar phenomena in different fields of science (5l-l7| . 

It seems to be the most challenging problem to treat 
Coulomb-like systems with high concentrations of inter¬ 
acting particles Q. When concentration increases, one 
can observe crystallization-like phenomena and transi¬ 
tion between different lattice symmetries [Ml. 

In current contribution we will treat concrete systems 
along with general formalism development. They were 
chosen as manifesting a lot of interesting effects. For ex¬ 
ample, dusty plasma may serve as perfect media for the 
experimental investigation of classical fluids and solids 
along with colloids The second system is “elec¬ 

trons on the liquid-helium surface” [231 - 1^ . It may un¬ 
dergo Wigner crystallization and transition be¬ 
tween triangular and square lattices [H, More¬ 

over, homogeneous distribution of electron density is not 
always stable, and there are critical parameters when 
spatial structures, especially periodic deformations and 
multi-electron dimples, are formed 


But our original goal is some general description of 
particles arranged in a Bravais lattice when interparticle 
potential may be “catastrophic”. It is rather difficult, 
since traditional methods of statistical mechanics can¬ 
not be applied to systems with Coulomb-like interactions 
because of free-energy divergences. Moreover, partition 
function can be exactly evaluated only for few model 
systems of interacting particles in the thermodynamical 
limit [37I - I40II . Thus, basing on [^ IdlMi^l . we will use al¬ 
ternative approach that employs Hubbard-Stratonovich 
representation of the partition function [d^. Besides, we 


use approach involving probability distribution function 
for the whole lattice as in [I[ , but construct it from single¬ 
particle distribution functions with a special procedure. 
Namely, we will expand and significantly supplement re¬ 
sults obtained in |35l | for two-dimensional systems. New 
results for three-dimensional case, that is far more com¬ 
plicated than two-dimensional, will be obtained. 

To be more rigorous: our goal is to find expression 
for free energy of particles arranged in a Bravais lattice 
even if inter-particle potential is “catastrophic”. Gen¬ 
eral equations will be applied to two concrete systems: 
three-dimensional — grains in dusty plasma, and two-di¬ 
mensional— electrons on the liquid-helium surface. In 
first case we’ll minimize this energy to find out which 
lattice should be observed experimentally. And for elec¬ 
trons localization distance will be found. Besides, this 
will show that presented method allows easily reduce di¬ 
mension number from three to two. 

Present article is organized as follows. It can be contin¬ 
gently divided into two parts: development of technique 
and its application. Sections U and m focus on general 
expression for free energy for any Bravais lattice and any 
potential. The former shows that not only classical, but 
lots of quantum systems may be treated this way. Be¬ 
sides, an expression for “entropy part” of free energy is 
calculated. And in section m interaction energy for pari- 
cles arranged in a Bravais lattice is calculated. 

Both this sections are equipped with examples (ira 
III AHII Cl) that show how to use developed formalism. 
This examples are rather useful, because they are later 
reused when treating physical systems. 

Concrete applications of developed formalizm start 
with section uni Here crystal of dust particles is con¬ 
sidered. In uni it is shown that hexagonal close packing 
seems to be the lattice we expect to be seen in exper¬ 
iment. This is in agreement with experiments Q and 
computer simulations 0- 

System of electrons on the liquid-helium surface is an¬ 
alyzed in II VI In llVI localization distance is calculated and 
compared to classical result [l^. For low temperatures 
T — >• 0 they coincide. 

For not to be confused with lots and lots of compli¬ 
cated mathematical computations we moved them to ap- 
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pendices IVII Thus article mostly contains some physical 
explanations and results. But using references to appen¬ 
dices, one may find here rigorous exposition of presented 
ideas as well. 


I. STATISTICAL DESCRIPTION OF 
INTERACTING PARTICLES SYSTEM 


First of all let’s outline class of systems we will be 
able to consider. Since everything is clear with classical 
systems, we will give heed to quantum ones. 

In this section we aim some general consideration of in¬ 
homogeneous system of interacting particles Nor 

system, neither interaction potential will be specified to 
keep generosity and extendability. As result, we will get 
free energy in some general form and ready for “apply¬ 
ing” translation symmetry (section [IT]) . Of course, along 
this section we will do some assumptions about system 
under consideration. Together they will form restrictions 
on quantum systems we can treat with proposed method. 

In current model macroscopic states of the system are 
described by occupation numbers. Besides, we suppose 
that it’s Hamiltonian has form: 

H + ^'^Vss'nsUs'. (1) 

s s,s' 

Here Eg is the additive part of particle energy (usually 
it is kinetic energy, but it can be particle’s energy in 
external field as well.), s indicates particle state, 14^' is 
interaction energy between particles in states s and s ', Ug 
is the occupation number of state s. We neglect quan¬ 
tum correlations, so only essentially classical systems are 
considered. 

Partition function for this kind of system will be: 

Z exp(-/3iJ), 

{Ws} 

where summation is performed over all possible states 
{us} of the system. 

Now we use some properties of Gaussian integrals over 
auxiliary fields, i.e. Hubbard-Stratonovich transforma¬ 
tion dilil: 


exp 



E 

.« s' 


^ss' 



OO 



— oo 



d 
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E 

.« s' 


^ggWsV’s 


with Dip = ris d^fs!\J det(27r/3a;ss')- Second-order de¬ 
pendence on occupation numbers is avoided by carrying 
them to introduced field and partition function can be 
written: 


Z = i Dp exp 


'^{ipg - Peg)ng - ^ E 


For canonical ensemble number of particles can be fixed 
by using Cauchy equation: 

that leads to A-particle partition function: 


Zn = ^ f j Dp exp 


1 


s,s' 


-(A + l)lne 


HE [^exp(f(/?s - /Ses)]”". 

« {"4 


Performing summation over occupation numbers: 

PF{p,i) = ]^'^V-J-pgPg,+ ( 2 ) 

s,s' 

+ <5^ In (1 - -f (Y-h 1) InC- 

S 

Here type of statistics is incorporated into 6 — it equals 
to -1-1 for Bose-Einstein, 0 for Maxwell-Boltzmann and 
— 1 for Fermi-Dirac statistics. 

Now we have expression for free energy in auxiliary 
field representation with ^ = exp(/3/r) for chemical activ¬ 
ity. Equation ([2]) contains all information about probable 
states of the system and corresponds to the sequence of 
equilibrium states with regard to their weights. We are 
interested in asymptotic value of partition function, but 
want to avoid using perturbation theory. Thus domain 
is extended to complex plane and saddle-point method is 
applied. 

Dominant contribution is made by states satisfying ex¬ 
trema condition: 


S/3F 

6p 

Variation of @ yields expression for saddle-point states: 


1 


= 0 ; 


E 


1 — 


■ = N + 1. 


(3a) 

(3b) 


Eor a certain state expression 

^g-PSa+iipa 

from (I3bl) can be treated as an average occupation num¬ 
ber. Now saddle-point states (that may be interpreted as 
thermodynamically stable distributions) can be obtained. 












3 


Equation (|3a|) contains inverse matrix, that is quite 
inconvenient for us, because its calculation is a chal¬ 
lenging mathematical problem even for simple poten¬ 
tials 113 ■ This problem can be avoided if we 

use ([4]) and perform inverse transformation 

ips = 

s' 

^ ^ss' fefs' = Vss' fsfs' ■ 

s,s' s,s' 

Than free energy can be written 

s (5) 

+ (fV+l)ln^(/). 

In terms of canonical ensemble we can write from (|4]) 
ln^(/) = ^ E + Es) + In/. - ln(l + <5/.)] , 

S 

Es=Y.Vss'fs'. 

s' 

and substitute last expression into free energy to 
rewrite its expression without chemical potential. 

F[f] = J2fses + lj2^ss'fsfs' + 

s s,s' 

+ 4 E - (/«+ 

^ s 

'-V-' 

Pent 

First two terms in (|6l) are kinetic and potential energy 
respectively, the last one is a contribution due to entropy 
Pent (should be equal to zero if T = 0). 

Now let’s do a quick test of obtained results. If a grand 
canonical ensemble with fixed chemical potential is con¬ 
sidered, we can get from (g]) 

= e/ 3 (E«-A*,) -S’ 

which is generalization of the well-known distribution 
with a chemical potential 


fj,s = fj.- Es. 

It is obvious that the “saddle point” and “mean field” ap¬ 
proximations are equivalents in this case. If an ideal gas 
{fig = Ai) is considered, we can obtain classical statistical 
distributions. 

Before winding up this section let’s review the results. 
We aimed to find some general expression for free energy 
in a mean-field approximation. As the one, equation ([6]) 
should be pointed out (and in some sense ([7]) that can be 
used to obtain particles distribution in grand canonical 
ensemble). Besides, we have got some restrictions on 
the class of systems we will be able to consider: special 
form of Hamiltonian o and negligible quantum cross¬ 
correlations. 


A. Example: two-dimensional system of localized 

Fermi particles— electrons on the liquid-helium 
surface 

Before we shift our attention to potential energy of 
particles arranged in a lattice, let’s do some example. 
With this subsection we aim two things: give some idea 
on how to use the above formalism and secondly— cal¬ 
culate contribution due to entropy for electrons on the 
liquid-helium surface (will be used in I®. 

Electrons on the liquid-helium surface have two degrees 
of freedom only [^, . Though we set <5 = — I (Fermi 

particles) and d = 2 (system is two-dimensional). Then 
equation ([5]) can be written 

f d^pd^r 1 

J {2-KhY 2m e/5(pV(2m)-M(r)) 1 + 

P ff d^pd^r d^p'd?r' 

^ 2 JJ (27r;i)2 (27rfi,)2 

^ _ V (l^-^'l) _ 

-p )) -p 

J (27r?i)2 1 + 1 ’’’ 

In -p 1 

g-^(p2/(2m)-/i(r )) _|_ I 




^g/3(p2/(2m)-/i(r)) 



Here integration sign f means integration over the whole 
phase space. 

Since last expression is quite cumbersome we should 
simplify it somehow. As first approximation we suppose 

£. = e(p) = p^l(2m). 

Actually, in the presence of an external field, dispersion 
relation should be a bit different. But we know, that 
system of Fermi particles, we will be interested in, is 
highly degenerated . So quadratic form for dispersion 
relation looks quite reasonable. 

Using last equation and introducing thermal length 


At = \/2'K‘^h?j3 I'm, 
we get from © following expression 


(9) 


pif) = 


d^p 


{2'KhY -p 1 


= Trln 

/\rp 


1^1 + 


( 10 ) 


which connects chemical potential /r('r) and particles 
density p{'r). 

Now we return to simplification of ®. Performing 
integration of ® over momentum and taking (nni and 
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do]) into account one may get (see appendix IVI All 
F[p] = - JJ (fr(fr'V{\r — f'\)p{r)p{r')+ 



One may notice that in case of Bose statistics we will 
loose the Bose-condensation effects due to this integra¬ 
tion. But we aim to use this equations for Fermi particles, 
so there is nothing to be worried about. 

We know that electrons on the liquid-helium surface 
exist in forms of fluid or Wigner crystal [2^-[3lj. We 
are interested in case when electrons are strongly local¬ 
ized, i.e. their coordinates are quite determined. But 
since temperature differs from zero, we expect that par¬ 
ticle’s position will fluctuate near its equilibrium loca¬ 
tion. Though we assume following distribution function 
Psp for single particle (somewhat analogous to form fac¬ 
tor in QFT) 


Psp(r) 


(27rs2)d/2 




( 12 ) 


where s is dispersion, or localization distance by physical 
meaning, and d = 2 because system is two-dimensiona0- 

Equation (HID is normal distribution which seems to 
be quite reasonable. At least this distribution is valid 
for ground state of quantum harmonic oscillator which is 
the simplest approximation for particle fluctuating near 
its equilibrium location. 

Last assumption we do is that “Gaussian” m is very 
“sharp” and every particle lays outside localization radius 
of any other particle. This means that s is much smaller 
then other characteristic distances in this system. So free 
energy per one particle can be written: 


-Fsp = ^ // r'\)Psp{r)p{r') + 




TT^S^ 


(13) 


87r2s2;3 6/3A2 


Here Fent is contribution to free energy by entropy 
and temperature. First summand in F^nt is obtained 
from direct integration of p1^ (see last term in (fTTI) '). To 
get second summand (penultimate term in (1111)1 we use 
“sharpness” of m and approximate 


Li2 


0 , r < s; 
7 r^/6, r > s. 


(14) 


With equations (|T^ and (ED we finish our considera¬ 
tion. We fulfilled our goal to show an example of using 
section’s U formalism. Moreover, now we have expression 
for entropy part Fent of free energy m, which we use 
later in IIVI 


II. PARTICLES ARRANGED IN A LATTICE 

In section H] we found some expressions for free energy 
of inhomogeneous system of interacting particles ([6D. In 
this section we will be interested in second term of m, 
i.e. potential energy. More accurately, we aim to find 
method of potential energy calculation for system of iden¬ 
tical (from inter-particle’s potential point of view) parti¬ 
cles arranged in a Bravais lattice. During this section we 
will show how to use translation symmetry to get in some 
cases more convenient expressions for potential energy. 

Last remark before we dive in. We aim to find appro¬ 
priate expressions for both two- and three-dimensional 
cases. But to avoid reduplication all along the section HIl 
only three-dimensional case will be treated as far more 
complicated. For two-dimensional case we will only pro¬ 
vide results and short description where difference comes 
from. 

Let V designates one subset of exact cover of (fig¬ 
ure [ID with similar domains containing strictly one par¬ 
ticle (relative position of the particle supposed to be the 
same in every domainal. Then free energy of the particle 
localized near Bravais lattice site can be written as (we 
neglect kinetic energy in ([HD since particles are localized 
and potential energy should be bigger) 

Fint= JJJ JJJ V {\f-f'\) p{r) p{r')d^f'(Iba) 

V R3 

Fs = JJJ JJJ V{\f-f'\)p{r)p{r')d^r'd^r,{15h) 

V V 

Fsp = Fint — Fg + Fent- (15c) 

Equation ED needs some explanations. Interaction 
energy is presented here with Pint — Fg . This “splitting” 
is just another way of saying we want to integrate over 

\ V. If one collects both expressions together and uses 
///r3\v = ///r 3 - fffv set regular integration 

over \ V for r'. 

In other words, since particle is localized near lattice 
site, we take some “part of space” V around it. Then 
by integration we calculate interaction energy between 
particle in V and all the rest particles in \ V. 

One may get a good intuitive idea about integration 
domain R^ \ V through discrete case. When we calculate 


^ We aim to reuse this and some of subsequent equations for 
systems with different number of dimensions. In such cases we 
will write them as functions of d (dimensionality) even if it is 
redundant at the point we introduce the equation first time. 


^ For example, Wigner-Seitz cells provide this kind of covering. 
But they often have too complicated shape and we will use par¬ 
allelepipeds instead (Figure 0. 
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Coulomb energy for i-th particle we write eiejjrij. 

Integration over \ V is somewhat analogous to j ^ i 
for discrete case. 

We will further call Fint (IlSal) interaction energy and 
Fs (I15bl) self-interaction energy. This “splitting” has 
one more convenience. Basing on physical background 
of considered potential we may want or not want to com¬ 
pensate self-interaction. Obviously, Coulomb potential 
should have self-interaction compensated — potential en¬ 
ergy is zero if there is only one particle in a system 
(considered later in example III Bll . On the other hand, 
self-interaction for effective potentials caused by influence 
through medium should be uncompensated, because even 
in absence of other particles medium will influence this 
one. Surface distortion is a good example. If particle has 
distorted a surface it has already done some contribution 
to system’s free energy independently on other particles 
(considered later in example III Cl) . Of course, it interacts 
with distortions created by other particles. But at the 
same time it interacts with “its own” surface distortion 
and this may be treated as “self-interaction”. 

Now lets consider (flSl) more carefully. All 3D Bra- 
vais lattices can be treated as constructed from paral¬ 
lelepipeds. We will use this property and use paral¬ 
lelepiped as integration region V. 



y 


Figure 1: One cell from exact cover V. Vectors a, b and 
c are base vectors for Bravais lattice. Angle between a 
and b is supposed to be a. Angle between c and XY 
plane is /3c and its projection on AF plane and a is ac- 
One particle is shown in the center. 

Next thing we want to do is expressing probability 
distribution function p for the whole through single¬ 
particle probability distribution psp. Basing on the in¬ 
verse lattice vectors fappendix IVI Bl) we will provide fol- 


lowing decomposition 

Pir) = Y. PkW)’ 

(16a) 


Pk = Pjjj fi(7')Psp(r )dV, 

(16b) 

V 

f,-(r) = 

(16c) 


/1 cot(a) 
a 

csc(a) 
b 

\0 0 



cot(/3c) sin(a — ac)\ 
a sin(a) 


cot(/3c) sin(ac) 

6 sin(a) 

csc{/3c) 

c 


(16d) 


/ 


Here * designates complex conjugation. It may be no¬ 
ticed that all vectors are treated as column vectors, p 
is mean particle density here. To get better intuitive 
understanding how to compute (HU one may check ex¬ 
ample [HX] for Psp from (fT^ . 

It may be shown that presented series are actually 
Fourier series for psp and p is periodic along vectors a, b 
and c 1appendix IVI Bl) . So we may treat p (r) as function 
on “composed” from single particle distributions “ar¬ 
ranged in a lattice”. This principle is demonstrated by 
Fig. [21 Since p has the same symmetry as the lattice and 




Figure 2: Two-dimensional analogy of transition from 
single-particle density probability function Psp{r) to 
many particle p(r ). Function p(r ) is defined for all R^ 
with two-dimensional analogy of (HU). 

locally is a good approximation for psp it will be treated 
as probability distribution function for the whole lattice. 

Matrix G comprises all geometry of the lattice under 
consideration. Performing minimization of Fgp over en¬ 
tries of G, one may find which lattice can be formed under 
certain conditions. We will use later this type of analysis 
for dust crystal in section Hill 

If we use equations m we can express (II5a|) in terms 
of Psp lappendix IVIDIl 

Fint = ^ IPfcl Vj:, (17a) 

^k=pJJJh ^ (17b) 

where p is mean particle density. 

Before moving further, one more thing should be em¬ 
phasized. This equations are useful because they allow 
us to deal with “catastrophic” potentials like Coulomb. 
One may notice that for any descending potential only 

CO 

^|popVb = 47rp J V{r)r‘^dr, 

0 

can be diverging and only if potential is “catastrophic”. 
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But this term does not contain “lattice geometry”, it 
depends on mean particle density only. This means, 
we can compare two lattices with equal mean particle 
density even if inter-particle potential is “catastrophic” 
(somewhat analogous to renormalization in QFT). 

With equations (II3, (HU and (USD we achieve our goal 
and finish consideration of interaction energy for parti¬ 
cles arranged in a lattice. Following examples show how 
to use obtained equations for screened Coulomb (exam¬ 
ple [110 and capillary interaction fexample III Cl) . They 
both use p obtained with pgp from (IT^ fexample lll All . 


liquid-helium surface. In both cases we have to deal with 
Coulomb potential — screened for dust particles and reg¬ 
ular for electrons. The latter will be treated as screened 
Coulomb with screening distance approaching infinity. 
During this example we aim to find Pint — Fs for screened 
Coulomb potential using methods of section |IT1 Consid¬ 
eration of physics of this two systems will be postponed 
until sections uni and El 
Screened Coulomb potential can be written as 


A. Example: calculating distribution function p for 
Gaussian single-particle distribution psp 


where Ad is Debye screening distance. It may be shown 
(appendix IVI Ell that expression for should be 


We suppose that all particles are arranged in a lat¬ 
tice and their positions are quite determined. But since 
temperature differs from zero we can expect that particle 
position fluctuates near its lattice site. Using considera¬ 
tions from oi we suppose it to be (fT^ with d = 3. 

Introduced constant s means dispersion of presented 
distribution. From physical point of view it can be 
treated as localization distance. Later, in section IIVI 
we’ll find its value and compare it with localization dis¬ 
tance obtained in different works. 

Using that Pspir) is a very “sharp” Gaussian, we ex¬ 
pand integration limits to in (I16bl) and perform inte¬ 
gration (appendix [VTCD 


As mentioned previously, we are interested in two- 
dimensional case as well. One may check, that there 
exists some sort of analogy of dSD for two dimensions 
[ 35 II . Difference will appeare in changing summations and 
integrations y IfI ri ^ f usage of dif¬ 

ferent matrix 


G = 


flja — cot(a)/a\ 


0 csc{a)/b 

V / 


(19) 


Performing calculations as in appendix IVI Cl with (fT^ 
and d = 2, one will get 


Pk = Pe 


-27i--’s^fc We k 


( 20 ) 


where p is two-, not three-dimensional density as in (1181) . 

Equations m and (HHD are aimed result of this exam¬ 
ple. Later they will be used in calculations of potential 
energy for Coulomb (IIIBI) and capillary (III Cl) interac¬ 
tions. 


B. Example: calculating Fi„t — Fa for screened 
Conlomb potential 




Airpq^ 


l/Xl + ATT^kTCG^k 


( 22 ) 


At this point we already have (l22l) and pj: (fTSl) . So 
we can write Pint from equation (I17al) 


F^nt = Y. 


-ATT^s^k GG^ k 




1/Xl+4Tr^ k^GG^k 


• (23) 


Compared to straight-forward method of computing, this 
one has better convergence. For big values of Ad we will 
see, that first summands in (1^ descend as l/|fcp. But 
if we write series containing interaction energy between 
this particle and any other, first summands will descend 
only as l/|fc|. 

In section [IT] we have mentioned that there is one com¬ 
pensatory term we need to take into account. This term 
will exclude Coulomb self-interaction. If we omit techni¬ 
cal details (see appendix IVI Fp explicit form of Fg can be 
obtained 


^ / '■''■'“•“’re')--'"*', (24) 

0 

or substituting screened Coulomb potential (|2T]) as V and 
approximating for the case when s is much smaller then 
screening distance Ad 


Fs = 



\/7rs Ad 


(25) 


Expressions we’ve got for Fi^t and Fg are a bit incon¬ 
venient. So we use that s is small compared to distances 
in the lattice (“sharp peak” approximation, section ED 
and perform some approximation. 

Eirst of all we introduce mean distance calculated from 
particle’s density 



In section Hill grains in dusty plasma will be considered 
and in section m we’ll set our eyes on electrons on the 


( 26 ) 
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and split sum (1231) into two parts. Further we will sup¬ 
pose lattice under consideration as not too much degen¬ 
erated. In other words, we will assume that it may be de¬ 
scribed as some “deformation” of cubic lattice (appendix 

ivrct 


Pint — Fg — 




kTG-^TG-lk/p^ 



PkTGGTj: 



-b 47rp(7^A|). 


(27) 


One may see that last term is equal to energy in case 
of evenly distributed charge. Other terms are corrections 
due to charge pointness and comprise all information on 
lattice geometry. 

It is very useful to notice, that this is the only sum¬ 
mand that approaches infinity if Xd —t oo- This means 
that (EZl) can be used for different lattice comparison even 
without screening. We will only need to measure energy 
starting from Airpq^X^, as “ground-level”. 

One more thing to notice. Coulomb is a long-range 
interaction and localization distance s should not play 
any significant role in approximation. Here we see that 
it is really absent. Rather different may be a picture 
for short-range interactions or effective interactions with 
uncompensated Fg ('example III Cl) . 


C. Example: calculating Pint for capillary 
interaction 


Last example we will consider is a rough estimation of 
capillary interaction energy. This results will be used 
in section IlYl More discussion on physics of the phe¬ 
nomenon and references we postpone until IlYl 

It is known, that electrons deform liquid-helium sur¬ 
face when pressed against with electric field and this 
causes effective interaction with potential [s^ 


V{r) = -^Koir/lo), (28) 

where qE is force caused by electric field, a — sur¬ 
face tension of liquid helium, r — distance between 
p articles, Kn— modified Bessel function, and Iq = 
y/a/gpHe 0 — capillary length that depends on the 
fluid properties only. 

Since we know that capillary interaction is essen¬ 
tially two-dimensional (particles should be localized on 
some surface, no three-dimensional analogue) we will 
use (HI) with ()19p and ()20l) . Performing calculations (ap¬ 
pendix |vIh|) we get 




q'^E'^ p 

^ 1/Iq k 


( 29 ) 


So we can immediately write 


Fint 




( 30 ) 


For our purposes (see section ITVT) a very rough approx¬ 
imation would be enough. So we use the same trick as in 
appendix IVI Gl but perform integration instead of sum¬ 
mation and dealing with theta functions (see appendix 

ETH) 


Fint = In {s/lo) ■ (31) 

dTTCT 

One should notice, that s <C ?0i so logarithm value is 
less then zero and thus whole expression m is less then 
zero. It coincides with the fact that capillary interaction 
makes particles attract each other [32 h34| . 


III. GRAINS IN DUSTY PLASMA, DUST 
CRYSTAL AND ITS LATTICE 

It is known that grains in dusty plasma interact and 
even show self-organization in form of melting and crys¬ 
tallization of dust crystalls [Ml. 

We aim to apply methods developed in section [IT] 
to this system. As a first approximation we suppose 
that grains interact only as charged particles (screened 
Coulomb potential). Obtained equation ([27|) (example 
III B|) allows to calculate energy of any lattice just insert¬ 
ing correct matrix G and performing summation. But 
much more interesting is finding the lattice with minimal 
energy when particles density p is fixed. Obtained result 
can be verified experimentally. 

Both sums in (l27l) are highly convergent and this is 
very convenient for numeric calculations. Thus we will 
limit summation with |fc| < 4. This limit was found 
experimentally: taking more summands does not change 
calculated minimal lattice and free energy. On the other 
hand, taking less results in wrong lattice. 

And one more thing to mension. Searching for lattice 
with minimal energy does not need calculation of full 
expression for — Eg. We need to minimize only parts 
with G. Moreover, they both contain g/Z, which means, 
that results of minimization do not depend of particle 
charge. 

Minimization was performed with standard function 
FindMinimum of Wolfram Mathematica v.9. We use (EZD 
as function to be minimized. Particles density is fixed in 
a following way. We consider 

P ^ o6csin(a) sin(/3c) ’ ^ 

suppose a = Z(l + da), b = 1{1 +6b) and express c through 
other parameters to keep constant charge density. Per¬ 
forming this one may notice that I cancel out in expo¬ 
nents and denominators (HZl). 
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Figure 3: Same lattice in different representations. It 
has translation symmetry with respect to two vectors a 
and b. Upper figure shows different choice for second 
translation vector b and b'. Using reflections one may 
show that lattice on bottom figure is the same as on the 
upper one. So b" is possible choice as well. When 
performing numerical calculations additional 
restrictions are used to overcome this uncertainty. 


Before considering results, one more remark should be 
done. Mapping between set of all possible parameter val¬ 
ues (translation vectors) and all possible Bravais lattices 
is not a bijection. It is rather a surjection— one lattice 
may be described with different sets of parameters (dif¬ 
ferent translation vectors). For minimization algorithms 
to be not confused, we provide some restrictions. 

Figure [3] shows uncertainty when choosing parameters 
for lattice description for 2D case. One may see that it 
can be eliminated if we claim that projection of & on a is 
less then a half of a. In terms of parameters this may be 
expressed as 


0 < (1 -I- db) cos(a) < (33) 

Restrictions for 3D case are more sophisticated. But 
they may be achieved similarly as (l33l) . Actually (l33l) 
should be kept for 3D as is. We only need to add analo¬ 
gous considerations for c. 

First of all we introduce projection of vector c on XY 
plane (see figure [J) and take particles density into ac¬ 
count (substituting c from (1331) 1 


Cx = ccos(/3c) cos(ac) = 
Cy = ccos(/3c) sin(ac) = 


I cot(/3c) cos(ac) 

(1 -I- (5a)(l -I- 6b) sin(Q;) ’ 
I cot(/3c) sin(ac) 

(1 -I- da)(l -I- Sb) sin(Q;) 


Obviously this projection should fall into parallelogram 
constructed on a and b. This condition is an obvious ex¬ 
tension of two-dimensional case to third dimension. But 


Figure 4: Projection of the same lattice in different 
representations on XY plane. Particle from layer above 
XY plane is shown with circle and below the plane with 
filled circle. Particles in XY plane are not explicitly 
shown. If particle from above the plane is projected into 
the upper half of ab parallelogram a reflection with 
respect to XY plane can be used to get representation 
where this particle is in lower half of parallelogram. 


we can strenghten this condition by using reflection with 
respect to XY plane (figure |4]). We can always assume 
that c is projected into lower part of parallelogram con¬ 
structed with a and b. 

As result we get following conditions 

^ Cy (1-I-(5&) sin(a) 

0 < -A < f- 1 -LA 

c„ cot(a) Cr Cy COt(Q;) _ ^ , 

^ I 

Now we can present results of numerical minimization 
in form of figure [S] and following table. 



Figure 5: Part of lattice obtained from minimization. 

As is seen, 3-D lattice consists of interleaving planes 
with 2-D hexagonal (triangular) lattices in each. The 
only difference between planes is shift: particles from 
one plane project onto “empty spaces” in another. 

Translation vectors are shown as well. 

Looking at figure [5l one my see, it is very similar to 
hexagonal close packing (HCP) lattice. HCP consists of 
two interleaving planes with triangular 2-D lattice in each 
as well. So we may check, if this is really HCP lattice by 
comparing their parameters. 
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Parameter 

Calculated 

value 

Hexagonal close 
packing (HCP) 

a 

ft! 1.04717 

7 r /3 ft 1.0472 

G.C 

ft 0.523589 

7 r /6 ft 0.5236 

/3c 

ft 0.955205 

arcsin ^2/3^ ft 0.9553 

a 

ft 1.122462 

^ ft 1.12246 

b 

ft 1.122462 

^ft 1.12246 


Results show that lattice obtained by numerical calcu¬ 
lations is really HCP. This coirmlies with experimental 
data and computer simulation . 


IV. ELECTRONS ON THE LIQUID-HELIUM 
SURFACE AND THEIR LOCALIZATION 
DISTANCE 


In this section we consider electrons on the liquid he¬ 
lium surface, calculate localization distance and compare 
it to obtained in other works. 

In the presence of an external field electrons can be 
pressed against the helium surface with significant force. 
But they cannot go throu gh t his surface, because quan¬ 
tum effects push them out . Thus they deform 

surface and this deformation, in its turn, changes interac¬ 
tion potential between electrons. The problem of finding 
explicit expression of this potential is solved by adding to 
Coulomb interaction effective capillary interaction. The 
lateral for two electrons on the helium surface was calcu¬ 
lated in [s^. 

We aim to find localization distance for electrons on the 
liquid helium surface. To achieve this goal we will mini¬ 
mize single-particle’s free energy Fgp over s (localization 
distance, see example llAI) . So first of all we construct Fgp 
and find dFsp/ds. 

In section U we developed a method for treating some 
quantum systems. Example II Al shows how to use it for 
two-dimensional system of Fermi-particles. So we con¬ 
sider Fent from m as free energy of one particle caused 
with entropy in the system. 

Potential energy is more complicated. In section |TT] 
we developed methods for its calculation. There are 
two parts of potential en ergy for electrons on the liq¬ 
uid helium surface ^t\. Coulomb interaction and 

capillary interaction. For Coulomb interaction we can 
find a shortcut — it was shown in [35l| that this part ap¬ 
proximately does not depend on localization distance s. 
Besides, one may think on approximate expression for 
Coulomb interaction from example III Bl (l27l) . It does not 
contain s, and the same behavior we may expect from 
Coulomb interaction for 2D lattice. Since we are inter¬ 
ested in derivative dFgp/ds, consideration of this part 
of potential energy may be omitted. More interesting is 
capillary interaction presented in example III Cl approx¬ 
imation m- It contains explicit dependence on s and 
thus should be taken into account. 


Bringing it all together we get 

dFgp TT^s q^E'^ 

ds 47r2/3s3 3/3A|, 4cr7rs’ 

independently on lattice we are considering. 

Substituting electron charge 9 —>■ e and solving 
tion dF/ds = 0 we get 

2 ( r Ibcr^TT'^ \ 

ATT^am y y j 

If temperature tends to be small T —^ 0 (e.g. /3 —>■ 


2TTah? 
m(?E'^ 

This result does not differ from the classical one pre¬ 
sented in p^ . 

If we consider (IMl) more carefully, some interesting re¬ 
sult may be obtained. Let’s consider low densities so that 
any value of s may be treated as much lower than inter¬ 
electron distance 1. Then we may notice that solution for 
s exists if 



equa- 

(34) 
00) 

(35) 


E'^ ^ Aair^k 


2100 


dyn 
cm ■ K 


(36) 


Last equation means that there is no lattice if tem¬ 
perature T is too high, or if field E is too weak. This 
statement is in agreement with our physical intuition. 


V. CONCLUSIONS 

In conclusion we want to make some overview of results 
we have obtained. 

In section [I] we started with a statistical description 
of inhomogeneous system of interacting particles. It was 
shown that variety of quantum systems may be described 
with no much difference from classical ones. We found 
an expression for free energy ©• The main feature of 
this equation appeares in expression for “entropy part” 
of free energy. Obtained result is applicable to a wide 
variety of systems with different number of dimensions. 

To make a step aside from pure theory we considered 
an example II Al of two-dimensional system of Fermi par¬ 
ticles. This result was later reused in section IlYl 

In section [Tl] we stated that particles are arranged in 
a lattice. It seems, that adding translation symmetry to 
consideration we can get some simplification for poten¬ 
tial energy expressions (uni and m- One may notice 
that for any “catastrophic” potential there is only one 
diverging term. But this term does not contain informa¬ 
tion on “lattice geometry”, it depends on mean particle 
density only. This means, we can compare two lattices 
with equal mean particle density even if inter-particle 
potential is “catastrophic”. 
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In this approach all potentials are classified into two 
groups: ones with self-interaction (often effective interac¬ 
tions through medium) and without. Examples III Bl and 
III Cl demonstrate how to use proposed technique in both 
cases. Later they are reused for treating grains in dusty 
plasma (section IIIII) and electrons on the liquid-helium 
surface 1 section HV)) . 

In section Iml we considered grains in dusty plasma. 
We started with potential energy obtained in examle fll Bl 
and minimized it over all possible lattices. Since there 
are no good approximations for theta-function it was too 
complicated to perform this minimization analytically. 
Thus we used numerical computations. This compu¬ 
tations have some features presented in a real-life case 
uni Results obtained in Ell show that we should expect 
hexagonal close packing (HCP) lattice for grains in dusty 
plasma. This coincides with numerical simulations in |10l | 
and experiments Q. 

Last section llVI shows how to use developed methods 
for electrons on the liquid-helium surface. It is shown 
that in some cases much more simple (but more rough 
as well) compared to section IIIII approximation may be 
obtained. This approximation is used to calculate lo¬ 
calization distance and then compared to classical result 
[ 13 . When temperature T —>• 0 they fully coincide. Be¬ 
sides, basing on this results some “guess” when lattice 
can exist is done. Obtained equation can be verified ex¬ 
perimentally. 


we would get 


mtA = 


m 


d‘‘fd"‘r'V{\f— r'|)x 


X In (1 -b In (l -b -b 


2 TTh?f 3 J 


d?r 


/3/r(r)ln(l + e^'"(”)) - 
- Li. ^ 


(37) 


which can be simplified further. 

First of all we will use Landen’s identity [45 


Li. (I - z) -b Li. ( 1 - 3 ) = -- In (z) 
z S C \ ]— 00 ; 0] 


1 


(38) 


with z = I -b e ^^0) amt expand In^ (l -b e ^^0)^, This 
leads to 


PF[d\ = 


d^rd^r'V{\r-r'\)x 


X In (l -b In (l -b -b 


27r/i^/3 , 


d^r 


Li. 


1 


1 -b 


TT 

Y 


+ iln2 (^I + 


VI. APPENDICES 


Now we substitute m and into the last expression 
and get (ITTl) . 


During this article we tried to provide rigorous expo¬ 
sition of presented ideas, but without physical sole being 
lost in lots of equations. Thus, calculations that heavily 
rely on mathematical transfomations only, were moved to 
appendices. Reader, interested only in “physical part”, 
may omit them. 


A. Integration over momentum in dF[fj,]. 

We start with ([5|). First two summands can be easily 
integrated if we mention 


1 


-b I I -b e-^ ’ 

but last summand should be integrated in terms of special 
functions, namely dilogarithm |45ll | 


Z 

U2iz) = - j 


ln(I - z) 


dz. 


If we mention that 


I -b 


= 1 - 


B. On the properties of fg(r'). 

In (fTBI) we present how to calculate probability distri¬ 
bution function for the whole system basing on single¬ 
particle probability distribution function. It is based on 
the decomposition of psp using fj:{r ) (libel) as basis func¬ 
tions in inverse lattice space. Here we add some con¬ 
sideration about connection between Fourier series and 
presented one. 

Suppose we have a function gif') defined on the unit 
cube U = [0; 1] x [0; 1] x [0;1]. This function can be 
expressed in terms of Fourier series 


9in = E 


feez3 


9k 


^(f)e-2^i(fc^r)^3- 


(39) 


1 -b e" 


Let’s consider a bijection G (llbdl) from V to U and its 
inverse G~^ (see Fig. [T] for geometrical reasoning) 

6 cos(a) ccos(ac) cos(/3c)^ 

0 6 sin(a) csin(ac) cos(/3c) j (40) 
yO 0 csin(/3c) 


G-i = 
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It may be seen that rows of G are components of in¬ 
verse lattice basis vectors. We are interested in case 
g{f) = Psp (^G~^r^. Since we know that Vr S V : 

Psp (^G~^Gr^ = Pspir) one can immediately write 

(41) 


And same way we consider second equation from (1391) 


9k = 


///..p(g-f) 


-2'Kik'^GG 


-.A 


G-V 


(42) 


where J is Jacobian. Here G is treated as new vari¬ 
able with domain V. 


Jacobian J 


G- 


= a6c sin(Q;) sin(/3c) = 1/p is ac¬ 


tually the volume of V or inverse mean particle density. 
Designating g^ as p-^ in (HD) and (|^ . changing vari¬ 
able in G“^r —7> r and designating exponent in (HD) 

as fg (it can be seen that exponent in (l42l) is ft, where 
* designates complex conjugation) we immediately get 
equations dUD. Moreover, obtained result means that all 
properties of Fourier series can be applied to decomposi¬ 
tion (fT6)l . 

Last thing to mention are periodical properties of pre¬ 
sented series. We are going to expand domain to so 
that p{r) will be defined everywhere in the space. So 
now we need to explore behavior of this function. From 
(I16dp one may see, if ^ G Z, m G Z and n G Z: 


G (r + la + mb + nc] = Gf + Icx + mCy + ncz, (43) 


where Cx, Cy and Cz are unit vectors along coordinate 
axis. With regard to (llBcIl 

(^r + la + mb + nc^ = f^(r ) (44) 

and from (I16al) 

p J- la + mb + nc'j = p{r). (45) 

Last equation justihes our view of connection between p 
and Psp as it is presented on Fig.jH 


Cartesian coordinates changing multiple integral to the 
product of integrals 


Pk 


(27rs2)3/2 






Last expression can be integrated if we use following 
relation [43| 


+ 00 



2 2_i_ 

—p X -^qx 


dx = 


— OO 


P 


5? (p2) > 0. 


As result we will get 


Pk= pH 
1=1 


Changing product of exponents to exponent of sum and 
mentioning that GG'^k we immediately 

get dUl). 


D. On the expression of Fi„f 

Lets start with (I15a|) . Inner integral has infinite bor¬ 
ders so we may rewrite this expression as follows 

Fint = JJJ JJJ V {\r'\) p{r) p{r + r')d^f'd^r. (46) 

v R3 


Since p (r) is real it can be replaced with complex conju¬ 
gate p* (r) without changing Pint- Now we can substitute 
p from (|16al) and mention that {r + f') = (r) fj: (r') 
(ll6cD . So if we designate V/ as in (I17bp it can be written 

r™. = - j: PfVi III k (n fj, (?) d”? 

fcez^ k'el? V 

Last equation can be simplified if we perform integra¬ 
tion. From appendix IVI Bl we expect orthogonality of 
functions. Let’s mention from » fs(r)fl,(r) = 


C. On the expression of pg for Gaussian psp 

We suppose that Psp{r) is either 0 everywhere in ]R^\V, 
or at least negligibly small. If so, we can change integra¬ 
tion limits in (jIBbll to infinite. 

At this point we are interested in specific form of psp 
thus it is explicitly substituted into (I16bl) . Be¬ 
sides we designate K = 2'Kk^G and rewrite expression in 


JJJ fki^)fi:,(.'r)d^r = abcsm{a)sm{/3c)Si:j:,. (47) 

V 

Here 5^j:, = h^,k'Jky,k'^dk^,k’^ — product of three Kro- 
necker’s delta. 

Mentioning that a6csin(a) sin(/lc) = 1/p, where p is 
mean particle density (one particle per V, see fig. [T] for 
geometrical reasoning) we get (I17al) . 
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E. Vg for screened Coulomb potential 


F. Calculation of Fs for screened Coulomb 
potential. 


Let’s consider equations (I16cl) and (I16d|) in spherical 


coordinates 




kx sin (etc 

, - a) 

ky sin(ac) 

a tan(/3c) 

sin(a) 

5tan(/3c) sin(a) 


ky 

cot(a)y 

5sin(a) 

'^X 1 5 

a J 


We will start with equation (IlSbl) and assume that 
Psp is a very “sharp” function 1 section ITAI) . With this as¬ 
sumption we can change integration limits to infinite and 
this, in turn, will allow us to perform variable exchange 
as we did in (l46t 

^^=III JJJ V{\r'\)Psp{r)psp{f+r')(ff'(fir. (48) 

B3 B3 


and rewrite Vg from (I17bll 


Vj: = pJ drViry j 
0 0 


dO sin(^) X 


X 


277 

J cos(e)+5l co&{ip-5ip^) sin(e)) 

0 

7(0) ' 


Integrating over (p we get 

I{e) = (27r(?^r sin(6<)) . 


First of all we perform some mathematical transfor¬ 
mations with (HU (in Cartesian coordinate system) and 
write the following 


Psp(r)psp(r+ r') 


_J_ -[r+r'/2]fis'^-r''^/{4s^) 


Now we can substitute this expression into (l48l) and per¬ 
form integration over r. Writing result in spherical co¬ 
ordinate system and integrating over angles we immedi¬ 
ately get (IMl) . 

Since we know explicit expression for V m we can 
substitute it into (IM)) and perform integration using def¬ 
inition of the complementary error function 


Then we substitute last expression into equation for Vj: 


tXJ TT 

V^ = 2ttp j V{ry j 
0 0 

X Jq ^27rggr sin(6*)^ sm{9)d9dr. 


We will consider screened Coulomb potential (EU, 
which means integration over r can be performed. One 
may use following relation 


J e-‘^^J^{Px)x''+^dr = 
0 


2a(2/3)'^r (i/-k|) 

0F(a2 + p2y+3/2 ’ 


5J(a) > |S(/3)|,K(z/) > -1 


erfc(a:) = ( e ^ dt. (49) 

Vtt j 

X 

As result we will get 

F, = (l - f^e*'/^-erfc(s/Ao)) . (50) 

Vtts \ \d j 

Since we may obviously suppose that s is very small com¬ 
pared to screening distance, then s ^ \d- So we can 
approximate last equation and as result we will get (1251) . 


G. Approximation of — Fa for screened 
Coulomb potential 


and combine it with equation for Vg to get 




2'Kpq^ (1/A_d — 27riggcos(d)) s\n{9)d9 

J / 2\ 3/2 ■ 

0 ( (1/Ad - 27ri5g cos(0)) ^ + (27rg^ sin(6»)^ 


Changing variable t = cos(0) and performing integration 
over t we will get 

^ _ ‘j'rrpq^ _ 

1/A|j -k 47r^gl + ' 

Last expression is equivalent to (1221) . One may check this 
by expanding g^ and g'^. 


We start with expression (12^ for Fint- But one may 
see that it contains expressions s^GG^. Components of 
matrix in (E5)) are proportional to different products 
of inverse distances in a lattice, e.g. 1/a^, l/(a&), and so 
on (see (jl6dl) b Since s <C a, s <C 6 and s <C c (assump¬ 


tion about “sharp Gaussian”) we expect 


s^GG^ 


< 1 . 


This means, that exponent “starts acting” only for terms 
with very large |A:|. In this section we will try to rewrite 
equations to get even better convergence then we have. 
First of all we rewrite expression for Fint as follows 


d^int 




27TS 


fcez3 
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One may easily integrate this expression and check that 
it is equal to (1^ . 

We notice, that variable s takes values much smaller 
then mean distance (E51) as well as much bigger. Thus 
we split this integral into two. First integration we will 
provide from 27rs to I and second from I to infinity. 

We assume that lattice is not too degenerated (sec¬ 
tion [110, which means that angles a and /3c should not 
significantly differ from 7r/2. So we may expect elements 
of s^GG^ to be less than 1 for first integral and greater 
than 1 for second. If so, then second integration may be 
performed and we end up with the sum that converges 
much better than (l23l) . 


Fin.t — 




''int — ^ 

TT 




ds-\- 




feez3 


e(0;s2GG^) 

^s^/\l-l^{l/[27v\of + k'^GG'^k) 

l/\l+ iTr^ k^GG^k 


From (llbdil one may check with Sylvester’s crite¬ 
rion [5l| that G (I16dp is a positive-definite matrix. Ob¬ 
viously, we expect and their product GG^ to be 
positive-definite matrices as well. Since s takes only pos¬ 
itive values s^GG^ is positive as well. This means high¬ 
lighted sum in the last equation is a Riemann theta func¬ 
tion [s^ at point z = 0, so we appropriately designate it 

with 0 ^0;s^GG^y 

Using modular transformation [s^ l 

7r^/2 

0(z;A) = ^=0(A-iz;A-i), (51) 

V det A 

where d is number of dimensions (3 in this case) we get 


0 


(O; 


s^GG^ = 


r3/2 


s3 det G 


0 (0;G-^^G-Vs2) . 


We used the fact that detG = detG'^ and that mul¬ 
tiplication by is equal to multiplication by diagonal 
matrix which has all elements equal to s^. Explicit ex¬ 
pression for G~^ may be used from (l40l) . Besides, we 
may notice that 


detG = 


1 

abc sin(a) sin(/3c) 


= P, 


(52) 


and rewrite expression for F^t as follows 


fcez^ 27I-S 


S2 


47rp(3'2 


fcgZ3 


I{k) 

gSVA|,-F(l/[27rAc]W|G’^fc|') 

l/Al, -f 47r2fc^GG^fc 


Now we can perform integration and find / (jtj . Since 

result is very complicated we first introduce two helper 
functions 


HL) = 




1 


erfc 


/„(L) = e“l‘^ '^l/(’^^^)erfc 


2\/^A_d 
G- ifc 


(2^Ab) 


-k 


ah 

27rAD 


and present result in terms of this functions 


J = 0^ = (j){2TTs) — 

j(fc^o) =V^(/_i(0-/-i(2^s))/ 


4G-ifc 


+ v^(/-si(0 - /+i(27rs)) / 


4G“W 


Using precise expression for Fg dsni) one may see 


Fg = 2y/TTq^e^ ^^^(P{2 tts). 


So subtracting Fg from Fint simply means neglection of 
this term. 

Our previous calculations do not rely on any approxi¬ 
mation, but expression for / ^ 0^ is very complicated, 

so we need to perform one. Before we do, one should 
show it is acceptable to approximate this sum. Subse¬ 
quent calculations are divided into two parts: uniform 
convergence proof and approximation itself. 

Proof of uniform convergence. We will use Cauchy 
criterion for sequence of functions fi{x) in domain E to 
achieve this goal: Ve > 0 37V Vm > n > TV : Vx G 

^ Further Y^kai? ^ S'® s 

sequence for convergence checking and entries of G~^ as 
well as s, I and Xd are treated as variables from domain of 
admissible parameters. We have already examined prop¬ 
erties of G~^ and will not get into details again. 

First of all we use s < I < Xd- This relation holds 
for physical system under consideration and we will use 
it in all subsequent calculations. This means fail) > 
fai2Trs) and thus all terms in sum will be positive and 
we can avoid using absolute value. The second thing that 
appears is following bounding 


0 < 


fcGZ3\0 fcGZ3\0 ^ 


Then we use relation erfc(x) < e~^^ for x > 0 [s^ and 
replace erfc in f±i{l)- Expanding squares in exponents 
and performing simplification we get 


E ^(«)< U 

feez3\o fcez3\o 


exp 


|G-ife|= F 

P 471-2 A^ 

G-ifc 






































14 


Besides < 4 was used. 

Now we can use known relation minpu^^i ||j 4 i ?||2 = 

VAniin, where Amin is minimal eigenvalue of A* A and 
* designates Hermitian adjoint [^. We know that 
det G~^ = 1/p > 0 and thus G~^ *G~^ has nonzero eigen¬ 
values. Taking the minimal one Amin and designating 

g = VAmin we claim G~^k > g 


Now we should change multidimensional summation 
to one-dimensional. Since all summands are positive we 
can rearrange them in any convenient order. Thus we will 
perform summation over “cube surfaces”. First of all let’s 
designate K„ = {k G \kx\ <nf\ \ky\ <nf\ \kz\ < n\. 
Then total sum can be expressed as 


+ CXD 

feez3\o 1 


E '(i). 

fe6K„\K„_i 


where number of summands in ]K„ \ K„_i can be easily 
calculated as (2n -I- 1)^ — (2n — 1)^ = 24n^ -|- 2. Minimal 
length for index vector in this set of indices is n. 
Collecting together previous calculations we get 


E 


(*) 


< 


+°° (24n^-|-2) exp ' 


r 


gn 


Further proof is supposed to be obvious and thus we avoid 
it. Proof of uniform convergence for first summand in 
Pint can be performed much more easily by presented 
scheme and thus is avoided as well. Instead one useful 
remark shoud be done. Series converges better when g/l 
is bigger. Taking physical meaning of this variables into 
account we get following statement: if we project vectors 
a, b and c on axis X, Y and Z respectively (see figure [T]) 
and divide it by mean distance in lattice we will 
get measure of how good current series converges, or by 
physical meaning— how close is current lattice to cubic. 

Approximation. We take into account s <C ? Xd, 
neglect all terms containing 1/Xd and s/Xu in erfc, and 
approximate erfc(x —>■ oo) ^ e~^ /{x^/n) [s^. Obviously 
fa{l) ^ /a(27rs) and thus all / with 27rs arguments are 
neglected. Besides, we neglect all small terms in expo¬ 
nent and denominators. To make this procedure more 


clear we state 


G-A 


/I > I/Xd- 


It means screening 


distance is big enough not to feel lattice deviation from 
cubic. This leads us to equations 




I 



1 

7’ 






To get the first one, identity erfc(x) = 1 — erf(x) and 
approximation erf(x —>■ 0 ) ^ 2x/^/tt were used. 

Now approximate equation can be written ( 1221 ). 


H. Vg for capillary interaction 

Let’s consider two-dimensional version of (II 6 cl) and 
(I16d|) in polar coordinates 

f^(f) = 


9k 


' i -2 

fV/y. 


\ 6 sin(a) 


ky cot(a) 


and rewrite Vg from (I17bll 


OC ZTZ 

V^ = pJ drV{r)r j . 


Integrating over p we get 

I = 27rJo (27r5gr) . 

Then we substitute last equality and expression for V (r) 
from (E51) into equation for Vg. It is known [7^, if a > 0 
and 6 > 0 


xJo(ax)Ko(6x)dx = 


(53) 


With this equation we easily obtain 14,,m for capillary 
interaction part 

^ (54, 


14 = -- 


(T 47r2g2 + l/Z2- 


Comparison with (1^ shows that they are equal. 


I. Approximation of Pint for capillary interaction 


We start with expression (1501) for Pint- First of all we 
rewrite it as follows 

Then we change summation to integration 
f 7^2 ■ There is a method of integration for this kind of 
functions [sij. As result we get 


Pint — 


pq^E'^ 


oSGll 


2a 


se 


p ,„= 

LG 


'■III 


se 


-bGiI 


ds 


2 s ^ det Q 


For two-dimension density it is true det Q — p. Besides, 
we perform integration over s and get 

where Fq is gamma-function. Since w know that s <C 
4 last equation can be approximated. We use known 
approximation [s^ ro(x —>■ 0 ) ~ —7 — ln(x), where 7 
is Euler-Mascheroni constant. We neglect this constant 
and exponent. The rest can be written as (EU. 
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